Gravitational evolution of a perturbed lattice and its fluid limit 



M. Joyce 

Laboratoire de Physique Nucleaire et de Hautes Energies, Universite de Paris VI, 
4, Place Jussieu, Tour 33 -RdC, 75252 Pans Cedex 05, France. 

B. Marcos 

Laboratoire de Physique Theorique, Universite de Paris XI, Bdtiment 210, 91405 Orsay, France. 



in 
o 
o 



in 

(N 

> 

CN 

o 
in 
o 

6 
to 

in 
c3 



A. Gabrielli 

SMC-INFM & ISC-CNR, Dipartimento di Fisica, 
Universita "La Sapienza" , P.le A. Moro 2, 1-00185 Rome, Italy. 

T. Baertschiger 

Dipartimento di Fisica, Universita "La Sapienza" , P.le A. Moro 2, 1-00185 Rome, Italy. 

F. Sylos Labini 

"E. Fermi" Center, Via Panisperna 89 A, Compendio del Viminale, 00184 - Rome, Italy. 

We apply a simple linearization, well known in solid state physics, to approximate the evolution at 
early times of cosmological iV-body simulations of gravity. In the limit that the initial perturbations, 
applied to an infinite perfect lattice, are at wavelengths much greater than the lattice spacing I the 
evolution is exactly that of a pressureless self-gravitating fluid treated in the analagous (Lagrangian) 
linearization, with the Zeldovich approximation as a sub-class of asymptotic solutions. Our less 
restricted approximation allows one to trace the evolution of the discrete distribution until the time 
when particles approach one another (i.e. "shell crossing"). We calculate modifications of the fluid 
evolution, explicitly dependent on / i.e. discreteness effects in the N body simulations. We note 
that these effects become increasingly important as the initial red-shift is increased at fixed I. The 
possible advantages of using a body centred cubic, rather than simple cubic, lattice are pointed out. 
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In current cosmological theories the physics of struc- 
ture formation in the universe reduces, over a large range 
of scales, to understanding the evolution of clustering un- 
der Newtonian gravity, with only a simple modification of 
the dynamical equations due to the expansion of the Uni- 
verse. The primary instrument for solving this problem 
is numerical iV-body simulation (NBS, see e.g.0). These 
simulations are most usually started from configurations 
which are simple cubic (sc) lattices perturbed in a man- 
ner prescribed by a theoretical cosmological model. In 
this letter we observe that, up to a change in sign in the 
force, the initial configuration is identical to the Coulomb 
lattice (or Wigner crystal) in solid state physics (see e.g. 
0), and we exploit this analogy to develop an approx- 
imation to the evolution of these simulations. We show 
that one obtains, for long wavelength perturbations, the 
evolution predicted by an analagous fluid description of 
the self-gravitating system, and in particular, as a spe- 
cial case, the Zeldovich approximation Further we 
can study precisely the deviations from this fluid-like be- 
haviour at shorter wavelengths arising from the discrete 
nature of the system. This analysis should be a useful 
step towards a precise quantitative understanding, which 
is currently lacking, of the role of discreteness in cosmo- 
logical NBS (see e.g. 0,H,E1)- O ne simple conclusion, 
for example, is that a body centred cubic (bcc) lattice 
may be a better choice of discretisation, as its spectrum 



has only growing modes with exponents bounded above 
by that of fluid linear theory. 

The equation of motion of particles moving under their 
mutual self-gravity is [l| 
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Here dots denote derivatives with respect to time t, Xj 
is the comoving position of the zth particle, of mass rrii, 
related to the physical coordinate by = a(t)xj, where 
a(t) is the scale factor of the background cosmology with 
Hubble constant H(t) — —. We treat a system of N 
point particles, of equal mass m, initially placed on a 
Bravais lattice, with periodic boundary conditions. Per- 
turbations from the Coulomb lattice are described simply 
by Eq. Q, with a(t) = 1 and Gm 2 — ► — e 2 (where e is 
the electronic charge). As written in Eq. Q the infinite 
sum giving the force on a particle is not explicitly well 
defined. It is calculated by solving the Poisson equation 
for the potential, with the mean mass density subtracted 
in the source term. In the cosmological case this is ap- 
propriate as the effect of the mean density is absorbed in 
the Hubble expansion; in the case of the Coulomb lattice 
it corresponds to the assumed presence of an oppositely 
charged neutralising background. 

We consider now perturbations about the perfect lat- 
tice. It is convenient to adopt the notation x,-(t) = 
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R + u(R, t) where R is the lattice vector of the ith parti- 
cle (which we can consider as its Lagrangian coordinate) , 
and u(R, t) is the displacement of the particle from R. 
Expanding to linear order in u(R, t) about the equilib- 
rium lattice configuration (in which the force on each 
particle is exactly zero), we obtain 

u(R, t) + 2#u(R, t) = 4V 2?(R - R')u(R', t) . (2) 

R' 

The matrix T> is known in solid state physics, for any in- 
teraction, as the dynamical matrix (see e.g. 2]). For 
gravity we have XV(R ^ 0) = Gm(|f - 3-%^) 
(where is the Kronecker delta), and P^„(0) = 

-E^o^(R) Hi- 

From the Bloch theorem for lattices it follows that V is 
diagonalised by plane waves in reciprocal space. Defining 
the Fourier transform by u(k, t) = Er e _ u(R, t) and 
its inverse as u(R, t) = -k Ek e ik u(k, t) (where the 
sum is over the first Brillouin zone), Eq. @ gives 

u(k, t) + 2H(t)u(k, t) = — ©(k)u(k, t) (3) 

where £>(k), the Fourier transform (FT) of £>(R), is a 
symmetric 3x3 matrix for each k. Diagonalising it one 
can determine, for each k, three orthonormal eigenvectors 
e„(k) and their eigenvalues w 2 (k) (n = 1,2,3), which 
obey Q the Kohn sum rule En w„(k) = — AnGpo, where 
po is the mean mass density. 

Given the initial displacements and velocities at a time 
t = tg, the dynamical evolution is then given as 

u(R, t) = ± E k E„[G(k, to) • e n (k)U n (k, t) 

+G(k,i ) • e„(k)K(k,t)]e„(k)e lk - R (4) 

where U n (k, t) and V n (k, t) are a set of linearly indepen- 
dent solutions of the mode equations 

f + 2H f = -^f (5) 
a- 5 

chosen so that U n (k, t ) = 1, ?7 n (k, to) = 0, V n (k, to) = 0, 
V n (k,t ) = 1. 

In Fig. n are shown the eigenvalues of the dynami- 
cal matrix, for gravity, for a 16 3 sc lattice, determined 
numerically by applying the linearization to a standard 
Ewald summation of the gravitational force (see e.g. 0). 
For convenience the eigenvalues have been normalized, 

with e n (k) = — ±™Qp , and they are plotted, as a func- 
tion of the modulus k = |k|, normalized to the Nyquist 
frequency = tt /I, where I is the lattice spacing. This 
diagonalisation can be performed rapidly even for the 
largest lattices used in current cosmological NBS, but 
the figure remains essentially unchanged except for an in- 
crease in the density of the eigenvalues. The lines in the 
figure connect the eigenvectors along some specific cho- 
sen directions, making the characteristic branch struc- 
ture of the eigenvectors evident. It can be shown 




FIG. 1: Eigenvalues e n (k) for a sc lattice. The lines connect 
eigenvectors with k in the specific directions indicated. Note 
that the two acoustic branches are degenerate in the [1,0,0] 
and [1,1,1] directions. 

that 2? M „(k — * 0) = —kpkvAnrGpo (where k — k/k), so 
the branch with eigenvalue tending to —AirGpo is lon- 
gitudinal (in this limit). In the Coulomb lattice this is 
the optical branch, describing oscillations for k — > with 
plasma frequency co^ — Aite 2 no/m (where no is the elec- 
tronic number density). There are then also two acoustic 
branches with eigenvalues tending to zero as k — > and 
which become purely transverse in this limit. A strik- 
ing feature of Fig|T] is that there are eigenvectors with 
e n (k) < 0, which correspond to unstable modes for the 
Coulomb system, with solutions to Eq. (JSJ U n (k,t) — 
cosh(|w„(k)|t) and V n (k,t) = (l/|w„(k)|) sinh(|w n (k)|t) 
(taking a — 1 and to = 0). Thus the sc Coulomb lattice 
is unstable to perturbations, which is not an unexpected 
result: the ground state of this classical system is known 
to be the bec lattice [S| , and these unstable modes in the 
sc lattice correspond to instabilities towards such lower 
energy configurations. For the case of gravity, in a static 
universe, these modes are sinusoidally oscillating, while 
the modes e n (k) > describe the expected exponential 
instabilities. Note further that, since the Kohn sum rule 
can be written E n e n(k) = 1, the appearance of modes 
with e„(k) > 1 is only possible when there are modes 
with e n (k) < 0. Without calculation we can thus con- 
clude that a bec lattice will have only unstable modes in 
the case of gravity, and that e n (k) < 1. We will return 
to this point below. 

The damping term coming from the expansion of the 
universe modifies these solutions to Eq. JSJ . For the case 
of an Einstein de Sitter (EdS, flat matter dominated) uni- 
verse, for which H 2 (t) = ^gpi and thus a = (t/t ) 2/3 , 
we find 

T7 ~ _ a+(k)(t/t )^^+a-(k)(t/t )-^ 

CCn(k) + On (k) 

a£(k) + a„ (k) 
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where a±(k) = \y\ + 24e„(k) ± 1]. Thus for e„(k) > 
there are, as in the static case, both a growing and 
a decaying solution. For e„(k) < the solutions are 
all power-law decaying. For e„(k) < — there is a 
weak remnant of the static universe oscillating behaviour: 
a^(k) are then complex, and it is simple to show that 
the mode functions are a product of a power law (t/to)~~s 
and a sinusoidal oscillation periodic in the logarithm of 
the evolution time ln(i/io)- 

Let us now consider the case that the initial fluctu- 
ations contain only modes s.t. kl <C 1. We have then 
simply for each k the longitudinal mode ei(k) = k, with 
ei(k) = 1, and two transverse modes with zero eigen- 
values. Using the corresponding mode functions from 
Eq. 10, and Eq. (0J, a simple calculation shows that 



u(Rt) =u x (R,*o) + U||(R,to) U^ + Ui)' 1 



-vii(R,to)io 
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-V ± (R,t )3t 
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where we have decomposed the particle displacements 
and peculiar velocities (v(R, t) = r s ; — ifr, — au(R, t)) 
into an irrotational (curl-free) part ay (R) = £k( a (R)' 
k)k e lkR , and a rotational part ai = a — aii . Us- 
ing the definition of the peculiar gravitational acceler- 
ation g(R, t) = f j — = a[ii + 2Hu], it is simple to 
show, using Eq. that g(R, t a ) = 47rGp U|| (R, t ) — 

^ttUii (R, to). Using this expression in Eq. Q, the diS- 
citQ II 

placement of each particle with respect to its initial po- 
sition (i.e. u(R, t) — u(R, to)) can be written solely in 
terms of the initial gravitational field g(R, t a ) and the 
components of the initial peculiar velocity, vj_(R, t ) and 
V|| (R, to). It is then easy to verify that the solution in 
Eq. J7J corresponds exactly to that derived in [9(, from 
a linearization of the Lagrangian equations for a self- 
gravitating fluid, for the displacements of fluid elements 
with respect to their Lagrangian coordinates As 
discussed in Q there are several limits of this expression 
which correspond to the so called Zeldovich approxima- 
tion (ZA), which assumes a decomposition of u(R, t) 
into a product of a function of time and a single vec- 
tor field defined at R. The most commonly used form 
of this approximation takes Uj_(R, to) = = vj_(R, t) 
and U||(R, i ) = f vy (R, to)to- This corresponds to set- 
ting the coefficients of all but the growing mode in 
Eq. (0) to zero i.e. it imposes directly the asymptotic 
behaviour of the general solution. We then have simply 
u(R, t) — |g(R, to)t^{t/to)^ which is precisely the solu- 
tion used standardly in setting up initial conditions for 
cosmological NBS (e.g. Q). 

This result provides a direct analytical derivation ex- 
plaining precisely the well documented success (see e.g. 
[l(j) of the ZA in describing the evolution of cosmological 
NBS, in particular in "truncated" forms of the approxi- 
mation in which initial short wavelength power is filtered 
The eigenvectors and the spectrum of eigenvalues 



contain, however, much more than this fluid limit. The 
expression Eq. Q gives an approximation to the full early 
time evolution of any perturbed lattice, treated as a full 
discrete A-body system. It therefore includes all modifi- 
cations of the theoretical fluid evolution in its regime of 
validity, which extends up to the time when particles ap- 
proach one another (i.e. up to close to "shell crossing"). 
We will report elsewhere detailed comparisons in numer- 
ical simulations of this approximation with the ZA and 
its improvements. In the rest of this letter we consider 
the quantification of the discreteness corrections to the 
pure fluid limit described by our approximation. 

Assuming still an EdS universe, and that the initial 
perturbations are set up in the standard manner us- 
ing the ZA, as described above, it follows directly from 
Eq. (3} that w M (k, t) = Y, v A^( k ' *)^( k > *o), where 
-V(M) = £j£4(i) + ^V„(t)](6n)/*(6n)^ (the k de- 
pendences on the right hand side are implicit). The full 
linearised evolution is encoded in this matrix, which can 
be calculated straighforwardly for any given lattice once 
the eigenvalues and eigenvectors have been found. One 
can then determine directly e.g. the power spectrum (PS) 
of the displacement fields <S M1/ (k, t) = u^(k, t)ii* v (k, t). 
Given S one can then calculate, by the method devel- 
oped in ^3 > the PS of the density field for the full point 
distribution. For small displacements (compared to I) 
and neglecting the terms describing the discreteness of 
the lattice, it is a good approximation to use the conti- 
nuity equation which gives Sp(k, t) w — ik • u(k, t), where 
Sp(k, t) is the FT of the density fluctuation field. It fol- 
lows that P(k,t) w A 2 P (k,t)P(k,to) where A P (k,t) = 
E MI A*wWM) and P(k,t) cx \Sp(k,t)\ 2 is the PS 
of the density fluctuations. It is simple to verify that in 
the fluid limit discussed above (kl — > 0) one obtains, as 
expected, A 2 P (k,t) = a 2 (t). 

In Fig. |21 is shown this amplification factor A 2 p (k, t), 
divided by a 2 . The scale factor chosen is a = 5, a value at 
which typical NBS reach shell crossing. Deviations from 
unity are a direct measure of the modification of the the- 
oretical evolution introduced by the discretisation. Note 
that Ap(k, a) is plotted as a function of k, each point 
corresponding to a different value of k. The fact that the 
evolution depends on the orientation of the vector k is a 
manifestation of the breaking of rotational invariance by 
the lattice discretisation. The three different symbols for 
the points correspond to three different intervals of the 
cosine of the minimum angle 9 between the vector k and 
one of the axes of the lattice. We thus see that the largest 
eigenvalues correspond to modes describing motion par- 
allel to one of the axes of the lattice. For a N 3 lattice 
and N even, for instance, the largest eigenvalue, with a 
growth law oc a 106 , is a longitudinal mode with k = fcjy 
and k parallel to the axes of the lattice, which describes 
the motion of pairs of adjacent infinite planes towards 
one another. Also shown in the figure is an average of 
A 2 P (k, a) over 25 bins of equal width in k, both for the 
I6 3 lattice from which the points have been calculated, 
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FIG. 2: Amplification function Af>(k, £) for the power spec- 
trum, divided by the fluid limit amplification (a 2 ), at a — 5, 
for a sc lattice. See text for details. 

and for a larger 64 3 lattice. 

We thus see that there are qualitatively two kinds 
of effects introduced by the discretisation: (i) an aver- 
age slowing down of the growth of the modes relative 
to the theoretical fluid evolution, and (ii) a pronounced 
anistropy in fc-space. There are notably a small fraction 
of modes (approximately 2.5 %) with growth exponents 
larger than in linear fluid theory (which, for sufficiently 
large a, will always dominate the evolution) . We can con- 
clude, however, as foreshadowed in the discussion above, 
that this evidently undesirable feature of the sc lattice 
discretisation can be circumvented by employing a bcc 
lattice. The known stability of this configuration of the 
Coulomb lattice Q implies that the fluid exponent is in 
this case an upper bound for all modes (and that there 
are no oscillating modes for the case of gravity). Further 
the bcc crystal is more isotropic (and indeed more com- 
pact than the sc lattice, and thus we would expect 
the effects of breaking of isotropy to be less pronounced. 
The average slowing down of the growth of the modes, 
by an amount which depends on the time and the di- 
mensionless product kl (at a = 5, as seen in Fig. [21 a 



10% effect at half the Nyquist frequency), on the other 
hand, would be expected to be a common feature of any 
discretisation (e.g. using "glassy" configurations or 
the discretisation developed in [Tql. 

This analysis does not, of course, allow one to conclude 
fully about the role of discreteness effects in the longer 
time evolution of such simulations, ft should provide, 
however, a first step to the quantitative understanding 
of these effects. One important conclusion which we can 
draw is the following: for a given physical wavelength, the 
discrepancy between the fluid and full evolution grows, 
up to shell crossing, with time. Thus, for a given phys- 
ical scale, discreteness effects increase when the starting 
time of the simulation is decreased. For the particular 
case of the sc lattice it is clear that artificial collapses of 
infinite planes along the axes of the lattice, which grow 
more rapidly than even long-wavelength fluid modes, will 
be increasingly priveleged as earlier starting times are 
taken. This implies that at least one of the conditions 
for keeping discreteness effects under control in an NBS 
will be a constraint on the initial time (i.e. that the 
starting red-shift be greater than some value, given the 
discretisation scale). We note that the initial red-shift is 
not a parameter considered in discussions of discreteness 
effects in NBS in the literature (e.g. |J,|3)- 

We can extend our treatment easily to incorporate a 
smoothing of the gravitational force up to a scale e. Here 
we have taken pure gravity (i.e e = 0) as in most cos- 
mological NBS e -C I, which gives negligible modification 
of our results. Just as in the analogous condensed mat- 
ter system, the method can also be extended to higher 
order, ft would be interesting in particular to map at 
higher order this description of the discrete system onto 
the corresponding order of fluid Lagrangian theory, which 
has been explored extensively in the cosmological liter- 
ature (see e.g. 16], and references therein). Further it 
should be possible to use the approach presented here to 
understand better the nature of existing approximations 
which go beyond the simple fluid limit, for example those 
involving pressure terms associated to velocity dispersion 
(see e.g. [l6|,[l3 and references therein). 
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